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A new approach for solving the Bethe ansatz (Gaudin-Richardson) equations of the standard 
pairing problem is established based on the Heine-Stieltjes correspondence. For k pairs of valence 
nucleons on n different single-particle levels, it is found that solutions of the Bethe ansatz equations 
can be obtained from one (k + 1) x (k + 1) and one (n — 1) x (k + 1) matrices, which are associated 
with the extended Heine-Stieltjes and Van Vleck polynomials, respectively. Since the coefficients in 
these polynomials are free from divergence with variations in contrast to the original Bethe ansatz 
equations, the approach thus provides with a new efficient and systematic way to solve the problem, 
which, by extension, can also be used to solve a large class of Gaudin-type quantum many-body 
problems and to establish a new efficient angular momentum projection method for multi-particle 
systems. 
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It is well known that the pairing force, similar to 
that in the Bardeen-Cooper-Schricffer (BCS) theory of 
superconductors [l[ , as one of main residual interactions 
introduced to the nuclear shell model, is key to manifest 
ground state properties and low energy spectroscopy of 
nuclei, such as binding energies, odd-even effects, single- 
particle occupancies, excitation spectra, electromagnetic 
transition rates, beta-decay probabilities, transfer reac- 
tion amplitudes, low-lying collective modes, level densi- 
ties, and moments of inertia, and so on Unlike 
electrons in solids, the drawbacks of the application of 
the BCS theory and its extensions to nuclei are notice- 
able due to the fact that the number of valence nucleons 
under the influence of the pairing force is too few to be 
treated by such particle-number nonconservation (quasi- 
particle) approximations 

Exact solutions to the standard pairing problem was 
first obtained by Richardson, now referred to as the 
Richardson- Gaudin method Recently, extensions 

to the Richardson-Gaudin theory have also been made by 
using the Bethe ansatz methodology The advan- 

tage of the Richardson-Gaudin solutions lies in the fact 
that the huge matrix in the Fock subspace is reduced to 
a set of equations, of which the number equals exactly to 
that of pairs of valence particles involved. However, less 
attention had been paid to the Richardson's solutions of 
the pairing problem in realistic calculations mainly be- 
cause the non-linear Bethe ansatz (Gaudin-Richardson) 
equations (BAEs) involved are very difficult to be solved 
numerically, especially for large size systems. Though 
there were a number of authors showing their efforts 
in designing algorithms for solutions with promising re- 
sults [13|]-[18|, obviously efficient procedure for solving 
the problem seems still unclear. Thus, a simple and clear 
approach to the problem is in demand. 

The Hamiltonian of the standard pairing model is 



given by 



H 



3=1 



-3 "J 



33' 



(1) 



where n is the total number of levels considered, G > 
is the overall pairing strength, {ej} are unequal single- 
particle energies, hj = J2 m a \m a jm is the number oper- 
ator for valence particles in the j-th level, and Sj~ = 

E m (-) J " m «L«] -m = (S+)*) are pair creation 

(annihilation) operators. Since the formalism for even- 
odd systems is similar, in the following, we only focus 
on the even-even seniority zero case. According to the 
Richardson-Gaudin method, fc-pair eigenstates of (1) can 
be written as 



\k;0 = S + (x 1 )S + (x 2 )---S + (x k )\0), 
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where |0) is the pairing vacuum state satisfying |0) = 
for all j, Xi (i — 1, 2, ■ • • , k) are spectral parameters to 
be determined. It can then be verified by using the cor- 
responding eigen-equation that (2) is the eigenstates of 
(1) only when the spectral parameters Xi (i = 1, 2, • • • , k) 
satisfy the following set of BAEs: 
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where pj = — (j + l/2)/2, with the corresponding eigen- 
encrgy given by E n , k = £) i=1 

Actually, as shown by Heine and Stieltjes, there is 
a one-to-one correspondence between every set of the 
Gaudin-Richardson type equations (BAEs) and a set of 
orthogonal polynomials, called by us the extended Heine- 
Stieltjes polynomials. Roots of these BAEs are zeros 



2 



of the polynomials, which can be interpreted as stable 
equilibrium positions in two dimensional complex plane 
for a set of free unit charges in an external electrostatic 
field The link between Richardson's BCS pair- 

ing model for nuclei and the corresponding electrostatic 
problem was thus established [2(J. According to Hcinc- 
Sticltjes correspondence, for nonzero pairing strcnght G, 
the polynomials y[x) with zeros corresponding to the so- 
lutions of Eq. ((3]) should satisfy the following second- 
order Fuchsian equation: 



A(x)y"(x) + B(x)y'(x) - V(x)y(x) = 0, 



(4) 



where A(x) = Oj=i ( x ~ ^ e i) * s a polynomial of degree n, 
B(x) is the polynomial with 



B(x)/A(x) = J2 
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and V(x) is called Van Vleck polynomials [19( of degree 
n — 1, which needs to be determined according to Eq. 

In the original electrostatic analogue considered by 
Heine and Stieltjes the parameters {pj} acting as 
fixed charges should all be positive with no external elec- 
trostatic field, 1/G — > 0. Therefore, the polynomials 
y(x) satisfying Eq. ((4]) with negative {pj} and 1/G ^ 
are thus called the extended Heine- Stieltjes polynomials, 
which tend to be the original Heine-Stieltjes polynomials 
with negative {pj} in the G — > oo limit. 

In search for polynomial solutions of (@J, we write 
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where {a^} and {bj} are the expansion coefficients to be 
determined. Substitution of © into Eq. ^ yields two 
matrix equations, the condition that coefficients in front 
of x 1 (i = 0, • • • , k) must be zero generates a (k + 1) x 
(k + 1) matrix F with Fv = b v, where the eigenvector 
v of F is just the expansion coefficients v = {oo, • ■ • , £Jfc}, 
while the condition that coefficients in front of x^ (i = 
k + 1, • ■ • ,n + k — 1) must be zero generates another (n — 
1) x (fc+1) upper-triangular matrix P with Pv = 0, which 
provides with unique solution of hi (i = 1, • • • , n — 1) 
in terms of {aj}. Entries of the two matrices are all 
linear with the coefficients {bi, 62, • • • , b n -\}. Matrices 
F and P can easily be constructed, for which a simple 
MATHEMATICA code is provided [21 1. 

Let the single-particle energies satisfy the interlac- 
ing condition e\ < ■ ■ ■ < e n . Real parts of zeros of 
y(x) satisfy the interlacing condition, —00 < Re(ii) < 
Re(x2) < • ■ • < Re(xfe) < +00, where Re(xj) lies in one 
of the n + 1 intervals (— 00, ei), (£1,62), (e 
and (e n ,+oo). It should be noted that many Re(xj) of 
adjacent zeros may lie within the same interval. When 



G —> 00, there will be only n intervals with (— 00, ei) 
being removed. The number of different such allowed 
configurations gives the possible solutions of y(x) and the 
corresponding V(x). The number of solutions of y(x), ex- 
cluding those with sum of zeros of y(x) complex, should 
equal to the number of levels produced by the standard 
pairing model, which is given by 
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where q = ~Y^ =1 p%- When pi = —1/2 for any i. which 
corresponds to the case of the Nilsson mean-field plus 
pairing model, rj(n,k) = n!/((n — k)\k\). Furthermore, 
if we set = 1 in y(x), the coefficient a^-i must equal 
to negative sum of zeros of y(x) with ak-i = —E n ^ = 
— X)*=i x i- Therefore, the solution corresponding to the 
largest real dk-i is that for the ground state of the sys- 
tem considered; those corresponding to the next largest 
real a^-i is that of the first excited state; and so on. In 
the standard pairing model, the solution with the same 
afc_i is unique except complex conjugation and permuta- 
tions within {xi}, which will be helpful in simplifying the 
calculation process, especially when only a few low-lying 
states are needed in the application. Since the coeffi- 
cients {a.,} and {bj} in F, P, and v are free from diver- 
gence with arbitrary variations in contrast to the original 
Bethe ansatz equations (3), one can use any standard re- 
cursive or iteration method to solve the problem with 
arbitrary initial values of these coefficients as desired. 
Because solving the eigen-equation Fv = 6ov, in which 
F is a (fc+ 1) x (k+ 1) matrix, is the only CPU time con- 
suming operation involved, the CPU time needed in the 
process should always be reasonable for k ~ 10 1 -10 3 and 
n ~ 10 1 -10 2 sufficient to realistic applications in nuclear 
physics. 

To demonstrate the new approach, we consider a 
simple example of k = 5 pairs in the sixth major shell 
with n = 5 levels, lh 7 / 2 , 2d 5/2 , 2d 3/2 , 3s 1/2 , and l/iu/ 2 , 
which is related to the application of the solution to Sm 
isotopes and difficult to be solved by directly us- 
ing the BAEs (3). We set single-particle energies to be 
equal spacing with = i, and the overall pairing strength 
G = 0.5. This example can now be dealt with easily by 
using the new polynomial approach, of which the number 
of solutions y(x) is 77(5, 5) = 71. First 5 sets of zeros of 
the corresponding polynomials y(x) and the correspond- 
ing coefficient ak-i are listed in Table I. The results can 
be obtained from the MATHEMATICA on a PC with 
in a minute. However, we observed there are some so- 
lutions with complex au-i, and some solutions are very 
close to each other like degenerate, which happened when 
MATHEMATICA built-in functions were used. Thus, 
the total number of solutions obtained is greater than 
77(71, fc), of which solutions with complex ak-i should be 
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discarded. Because there is little to be known about 
these polynomials with negative charges, further study 
is needed to see whether there are indeed solutions with 
complex asfc_i. Both complex a^-i and near degenerate 
issues may all be due to the Newtonian iteration adopted 
in the original MATHEMATICA package for solving a 
set of equations. Special codes designed suitable for the 
approach seem needed to overcome the ambiguity. 



TABLE I: First 5 sets of zeros of the possible solutions of 
the extended Heine-Stieltjes polynomials and the correspond- 
ing eigen- energies (in arbitrary unit) of the standard pair- 
ing model Hamiltonian (1) in the case of k — 5 pairs in the 
sixth major shell with n = 5 single-particle levels 1/17/2, 2ds/2, 
2d 3 / 2 , 3si/2, and l/in/2- The single-particle energies used are 
ti = i, and the overall pairing strength G = 0.5. 



Zeros of the polynomials 



Xl 


= -1.4993, x 2 = -1.1412 - 2.1396J, 


-3.6158 


X3 


= -1.1412 + 2.13961, x 4 


= 0.0829 - 4.5018i, 




X5 


= 0.0829 + 4.5018i 






Xl 


= -0.5078 - 1.0411i, x 2 


= -0.5078 + 1.0411i, 


3.0299 


X3 


= 0.5469 - 3.3066i, x 4 = 


0.5469 + 3.3066i, 




XS 


= 2.9517 






Xl 


= -0.9234 - 1.0718i, x 2 


= -0.9234 + 1.0718i, 


3.5444 


X3 


= 0.0573 - 3.3613i, x 4 = 


0.0573 + 3.3613i, 




XS 


= 5.2767 






Xl 


= -1.1244 - 1.0987i, x 2 


= -1.1244 + 1.0987i, 


4.8379 


X3 


= -0.1739 - 3.4422i, x 4 


= -0.1739 + 3.4422i, 




XS 


= 7.4346 






Xl 


= -1.2032 - 1.1109i, x 2 


= -1.2032 + 1.1109i, 


5.77020 


X3 


= -0.2619 - 3.4804i, x 4 


= -0.2619 + 3.4804i, 




X5 


= 8.7004 







In addition, as shown in our previous study [22J, 
a new angular momentum projection method for multi- 
particle systems can be established based on BAEs sim- 
ilar to (3). In fact, for n angular momenta ji (i = 
1,2, ■ ■ ■ ,n), the mult i- particle state with total angular 
momentum J = ^ i ji — k can be written as 

\ V , J, M = J) = J-(x x )J-(x 2 )---J-(x k )\h.w.), (8) 

where rj is a quantum number needed to resolve the 
multi-occurrence of J, h.w.) is the highest weight single- 
particle product state with mi = ji, ■ ■ ■ ,j n , m n = 
jn), and 



(9) 



and Ci (i = 1, 2, ■ • ■ , n) can be taken as any set of un- 
equal numbers 22| . Acting the total angular momentum 
raising operator J + to (JS} with J + \ri, J, M = J) = 0, 
one obtains the BAEs of {xi} the same as Eq. © in 
the replacement pi = — (ji + l/2)/2 in the G — > 00 limit. 
Therefore, once the solutions of (3) in the G — > 00 limit 
are obtained, the resultant {xi} once and for all deter- 
mine the mult i- particle state with good angular momen- 
tum J. The number of solutions of (3) equals exactly to 
the number of occurrence of J for the given system. The 
polynomial solutions (4) in this case exactly become the 
Heine-Stieltjes polynomials mentioned previously. This 
angular momentum projection is certainly much simpler 
than the projection operator technique Q and that based 
on the permutation group method [23j. 

In summary, we have established a new approach for 
solving the standard pairing problem based on a sound 
mathematical foundation — the extended Heine-Stieltjes 
polynomials and the corresponding Van Vleck polyno- 
mials satisfying the polynomial solutions of the second 
order Fuchsian equation. Thus, we reach the go of the 
Richardson-Gaudin theory via the Heine-Stieltjes corre- 
spondence, from which the exact solutions to the problem 
can practically be realized based on two matrix equa- 
tions. The approach can easily be extended and applied 
to solve a large class of Gaudin-type quantum many-body 
problems. A new efficient angular momentum projection 
method for multi-particle systems is thus proposed as a 
byproduct, of which the application to either boson or 
fcrmion systems will be studied elsewhere. 
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in which J. is the angular momentum lowering opera- 
tor only acting on the i-th single-particle state \ji, mj), 
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